!------------------------------------------------------------------------
!stomatal resistance (LSM_Bonan) added as optional calculation by AYu
!------------------------------------------------------------------------
      subroutine stomata (apar  ,&
     &                    ei   ,ea   ,tpbl   ,pgr  ,&
     &                    Z,t,wll, WSSL, WSSG,& 
     &                    PU, PV, PCH, rs, psn, psn_un)

!------------------------ code history ---------------------------
!source file:       stomata.F
!purpose:           leaf stomatal resistance and leaf photosynthesis
!date last revised: March 1996 - lsm version 1
!author:            Gordon Bonan
!standardized:      J. Truesdale, Feb. 1996
!reviewed:          G. Bonan, Feb. 1996
!-----------------------------------------------------------------
!------------------------ input/output variables -----------------
!input

      use INMSOIL_NUMPARAMS, only : &
     & ML, MS

      INCLUDE 'soil.inc'

!      PARAMETER (nv2=7) !number of vegetation types
      !USE PTRMODPHY, ONLY : NV2

      REAL (KIND=8)  ei          !vapor pressure inside leaf (sat vapor press at tv) (pa)
      REAL (KIND=8)  ea          !vapor pressure of canopy air (pa)
      REAL (KIND=8)  apar(nv2)   !par absorbed per unit lai (w/m**2)
      REAL (KIND=8)  pgr         !air pressure at agcm reference height (pa)
      REAL (KIND=8)  tpbl    !air temperature at agcm reference height (kelvin)
          REAL (KIND=4)  Z(ML)           ! soil layer depth,cm
          REAL (KIND=4)  T(ML)  !soil temperature degC
          REAL  (KIND=4) wll(ML) !water content (g/g)
          REAL (KIND=8) WSSL(ML,nv2) !lower limit water content for transpiration 
      REAL (KIND=8) WSSG(ML,nv2) !optimum water content for transpiration
      REAL (KIND=8) PU, PV !U and V wind speed (m s-1)
          REAL (KIND=8) PCH ! turbulent drag coefficient
!output
      REAL (KIND=8)  rs(nv2)     !leaf stomatal resistance (s/m)
          REAL (KIND=8) psn(nv2) !foliage photosynthesis (umol co2 /m**2/ s) [always +]
          REAL (KIND=8) psn_un(nv2) !non water-limeted foliage photosynthesis (umol co2 /m**2/ s) [always +]
!-----------------------------------------------------------------

!------------------------ local variables ------------------------
          REAL (KIND=8)   rb(nv2)     !boundary layer resistance (s/m)
      REAL (KIND=8)    btran(nv2)  !soil water transpiration factor (0 to 1)
          integer i,k,ik   !array/loop index
      integer iter     !iteration index
      integer niter    !number of iterations
          REAL (KIND=8)    tfrz        !freezing point (kelvin)
      REAL (KIND=8)    mpe         !prevents division by zero errors
          REAL (KIND=8)    o2               !atmospheric o2 concentration (pa)
      REAL (KIND=8)    co2               !atmospheric co2 concentration (pa)
          REAL (KIND=8)    igs(nv2)   !growing season index (0=off, 1=on)
          REAL (KIND=8)    tv(nv2)     !foliage temperature (kelvin)

     REAL (KIND=8)    qe25(nv2)   !quantum efficiency at 25c (umol co2 / umol photon)
     REAL (KIND=8)    ko25(nv2)   !o2 michaelis-menten constant at 25c (pa)
     REAL (KIND=8)    kc25(nv2)   !co2 michaelis-menten constant at 25c (pa)
     REAL (KIND=8)   vcmx25(nv2) !amax1imum rate of carboxylation at 25c (umol co2/m**2/s)
         REAL (KIND=8)  vcmx_un(nv2)   !non water-limited vcmax
     REAL (KIND=8)   aqe(nv2)    !q10 for qe25
     REAL (KIND=8)   ako(nv2)    !q10 for ko25
     REAL (KIND=8)    akc(nv2)    !q10 for kc25
      REAL (KIND=8)    avcmx(nv2)  !q10 for vcmx25
      REAL (KIND=8)   bp(nv2)     !minimum leaf conductance (umol/m**2/s)
      REAL (KIND=8)   mp(nv2)     !slope for conductance-to-photosynthesis relationship 
      REAL (KIND=8)   folnmx(nv2) !foliage nitrogen concentration when f(n)=1 (%)
       REAL (KIND=8)   foln(nv2)  !foliage nitrogen concentration (%)
      REAL (KIND=8)   c3psn(nv2)  !photosynthetic pathway: 0. = c4, 1. = c3
          REAL (KIND=8)    rdp(nv2)
          REAL (KIND=8) r_a(nv2), r_b(nv2) ! root functions
      REAL (KIND=8) prevent 



      data niter /3/
      save niter

!      REAL (KIND=8)  ab          !used in statement functions
!      REAL (KIND=8)  bc          !used in statement functions
      REAL (KIND=8)  fsto1          !generic temperature response (statement function)
      REAL (KIND=8)  fsto2         !generic temperature inhibition (statement function)
      REAL (KIND=8)  tc          !foliage temperature (degree celsius)
      REAL (KIND=8)  cs          !co2 concentration at leaf surface (pa)
      REAL (KIND=8)  kc          !co2 michaelis-menten constant (pa)
      REAL (KIND=8)  ko          !o2 michaelis-menten constant (pa)
      REAL (KIND=8)  a,b,c,q     !intermediate calculations for rs
      REAL (KIND=8)  r1,r2       !roots for rs
      REAL (KIND=8)  fnf         !foliage nitrogen adjustment factor (0 to 1)
      REAL (KIND=8)  ppf         !absorb photosynthetic photon flux (umol photons/m**2/s)
      REAL (KIND=8)  wc, wc_un          !rubisco limited photosynthesis (umol co2/m**2/s)
      REAL (KIND=8)  wj          !light limited photosynthesis (umol co2/m**2/s)
      REAL (KIND=8)  we, we_un         !export limited photosynthesis (umol co2/m**2/s)
      REAL (KIND=8)  cp(nv2)     !co2 compensation point (pa)
      REAL (KIND=8)  ci(nv2)     !internal co2 (pa)
      REAL (KIND=8)  awc(nv2)    !intermediate calcuation for wc
      REAL (KIND=8)  vcmx(nv2)   !amax1imum rate of carboxylation (umol co2/m**2/s)
      REAL (KIND=8)  ej(nv2)     !electron transport (umol co2/m**2/s)
      REAL (KIND=8)  cea(nv2)    !constrain ea or else model blows up
      REAL (KIND=8)  cf(nv2)     !s m**2/umol -> s/m
      REAL (KIND=8)  htop(nv2)    !top of the canopy (m)
      REAL (KIND=8)  dleaf(nv2), d(nv2), z0m(nv2)
      REAL (KIND=8) rl, udyn, uztop(nv2) 

!        REAL (KIND=8)  ROOTSM(nv2)
        REAL (KIND=8)  tmin(nv2)
        REAL (KIND=8)  r(ML)          !REAL (KIND=8) tive root fraction in soil layers
    REAL (KIND=8)  ww(ML)          !water content modifyer for evaporation
! minimum temperature for photosynthesis (kelvin)
      data (tmin(i), i=1,nv2) /278.16,273.16,268.16,273.16,&
     &           273.16,273.16,273.16/ !,  0.00/

!photosynthetic pathway: c3 = 1, c4 = 0
      data c3psn /7*1./ !, 1./

! co2 michaelis-menten constant at 25c (pa)
      data kc25 /7*30./

!q10 for kc25
      data akc /7*2.1/

! o2 michaelis-menten constant at 25c (pa)
      data ko25 /7*30000./

! q10 for ko25
      data ako /7*1.2/

! maximum rate of carboxylation at 25c (umol co2/m**2/s)
      data vcmx25 /50.,33.,33.,33.,17.,&
     &             33.,50./ !, 0./
!       data vcmx25 /7*250./
! q10 for vcmx25
      data avcmx /7*2.4/

! minimum leaf conductance (umol/m**2/s)
      data bp /7*8210./ ! 7*2000./

! slope for conductance-to-photosynthesis relationship
       data mp /9.,9.,6.,9.,9.,9.,9./ !,9./
!      data mp /7*9.0/ 
! quantum efficiency at 25c (umol co2 / umol photon)
      data qe25 /7*0.06/ !, 0.00/

! q10 for qe25
      data aqe /7*1.0/
!foliage nitrogen concentration when f(n)=1 (%)
      data folnmx /7*1.5/ !, 0.0/
! empirical parameter that defines root fraction decrease with depth
      data rdp /0.94,0.94,0.94,0.97,&
     &          0.97,0.94,0.94/ !,1.00/
!top of canopy (m)
      data (htop(i), i=1,nv2) /35.0,20.0,17.0, 0.5,&  
     &           0.5, 0.5, 0.5/ !, 0.0/

      data dleaf /7*0.04/
      data z0m /2.62, 1.1, 0.94,0.06, &
          & 0.06,0.06,0.06/
          data d /23.45, 13.4, 11.39, 0.34, &
          & 0.34, 0.34, 0.34/
          data r_a /7.0, 6.0, 7.0, 11.0, 7.0, 7.0, 6.0 /
          data r_b /1.0, 2.0, 2.0, 2.0, 1.5, 1.5, 3.0 /

          prevent=0.000001     
!--------------------------------------------------------------------


        
        o2=20950.0
        co2 = 38.308

        mpe = 0.000001
        tfrz = 273.15
        do k = 1,nv2
                tv(k) = tpbl
        end do
!-----------------------------------------------------------------
  

!-----------------------------------------------------------------
! Phenology: define growing season
! -----------------------------------------------------------------
      do k = 1, nv2
         if (tv(k) .gt. tmin(k)) then
            igs(k) = 1
         else
            igs(k) = 0
         end if
      end do
! -----------------------------------------------------------------

!-------------------------------------------------------------------
!set N to non-limited value 
!--------------------------------------------------------------------
       do k = 1, nv2
       foln(k)=folnmx(k)
           end do

! beta factor 'btran' calculation

    do k = 1,nv2
    btran(k) = 0.0
        do i=MS+1,ML-1
                   r(i)=0.5*(exp(-r_a(k)*z(i-1))+exp(-r_b(k)*z(i-1))- &
         & exp(-r_a(k)*z(i))-exp(-r_b(k)*z(i)))
               if (t(i)+273.15.gt.tfrz)then
               ww(i)=MIN(1._8,MAX(mpe,(wll(i)-WSSL(i,k))/(WSSG(i,k)-WSSL(i,k))))
           else
               ww(i)=0.01
                end if
        btran(k)=btran(k)+r(i)*ww(i)
        end do
        end do
        
!rb calculations (Bonan, LSM)
    rl=1./dmax1((PCH*sqrt(PU**2+PV**2)),prevent)
        udyn=sqrt(sqrt((PU/rl)**2+(PV/rl)**2))
!       write(*,*) udyn
    do k = 1,nv2
        uztop(k)=udyn/0.4*log((htop(k)-d(k))/z0m(k)) 
    rb(k)=1./dmax1((0.02/3.0*sqrt(uztop(k)/dleaf(k))*(1-exp(-3./2.))),prevent)
!         write(*,*) rb(k)
        end do
!initialize rs=rsamax1 and psn=0 because will only do calculations
!for apar > 0, in which case rs <= rsamax1 and psn >= 0

      do k = 1, nv2
         cf(k) = pgr/(8.314*tpbl)*1.E+6 
         rs(k) =1./bp(k) !5000._8/cf(k)
         rb(k)=rb(k)/cf(k)
!rb: s/m -> s m**2 / umol 
         psn(k) = 0.
      end do
  

!ci iteration

      do iter = 1, niter
         do k = 1, nv2
                 if (apar(k).gt.0.0)then
         fnf = dmin1( foln(k)/dmax1(mpe,folnmx(k)), 1.0 )
         tc = tv(k)-tfrz                            
         ppf = 4.6*apar(k)                  
         ej(k) = ppf*qe25(k)
         kc = kc25(k) * fsto1(akc(k),tc)       
         ko = ko25(k) * fsto1(ako(k),tc)
         awc(k) = kc * (1.+o2/ko)
         cp(k) = 0.5*kc/ko*o2*0.21
         vcmx(k) = vcmx25(k) * fsto1(avcmx(k),tc) / fsto2(tc) * fnf * btran(k)
         vcmx_un(k)=vcmx25(k) * fsto1(avcmx(k),tc) / fsto2(tc) * fnf  !AYu
!first guess ci
          
         if(iter.EQ.1) ci(k) = 0.7*co2*c3psn(k) + 0.4*co2*(1.-c3psn(k))  

!constrain ea

         cea(k) = dmax1(0.25*ei*c3psn(k)+0.40*ei*(1.-c3psn(k)), &
     &                  dmin1(ea,ei) )

            wj = dmax1(ci(k)-cp(k),0.)*ej(k)/(ci(k)+2.*cp(k))*c3psn(k) +&
     &           ej(k)*(1.-c3psn(k))
            wc = dmax1(ci(k)-cp(k),0.)*vcmx(k)/(ci(k)+awc(k))*c3psn(k) +&
     &           vcmx(k)*(1.-c3psn(k))
            we = 0.5*vcmx(k)*c3psn(k) +&
     &           4000.*vcmx(k)*ci(k)/pgr*(1.-c3psn(k))
            wc_un = dmax1(ci(k)-cp(k),0.)*vcmx_un(k)/(ci(k)+awc(k))* &
     &          c3psn(k)+vcmx_un(k)*(1.-c3psn(k)) !AYu
            we_un = 0.5*vcmx_un(k)*c3psn(k) + &
     &           4000.*vcmx_un(k)*ci(k)/pgr*(1.-c3psn(k)) !AYu

            psn(k) = dmin1(wj,wc,we) * igs(k)
            psn_un(k) = dmin1(wj,wc_un,we_un) * igs(k) !AYu 
               
            cs = dmax1( co2-1.37*rb(k)*pgr*psn(k), mpe )
            a = mp(k)*psn(k)*pgr*cea(k) / (cs*ei) + bp(k)
            b = ( mp(k)*psn(k)*pgr/cs + bp(k) ) * rb(k) - 1.
            c = -rb(k)
            if (b .ge. 0.) then
               q = -0.5*( b + sqrt(b*b-4.*a*c) )
            else
               q = -0.5*( b - sqrt(b*b-4.*a*c) )
            end if
            r1 = q/a
            r2 = c/q
            rs(k) = dmax1(r1,r2)
!            write(*,*) iter,igs(k),apar(k), psn(k), rs(k)
            ci(k) = dmax1( cs-psn(k)*pgr*1.65*rs(k), 0. )
            end if
            end do
      end do

!rs, rb:  s m**2 / umol -> s/m 

      do k = 1, nv2
       
         rs(k) =rs(k)*cf(k)
         rb(k) = rb(k)*cf(k)
          
      end do

      return
      end

!--------------------------------------------------------------------------
      function fsto1(ab, bc)
          REAL (KIND=8) fsto1
          REAL (KIND=8)  ab, bc
      fsto1 = ab**((bc-25.)/10.)
          end function fsto1

          function fsto2(ab)
          REAL (KIND=8) fsto2
          REAL (KIND=8)  ab      
          fsto2 = 1. + exp((-2.2e05+710.*(ab+273.16))/(8.314*(ab+273.16)))
          end function fsto2
 
